\section{MCMC}

{\color{red}\textbf{[历史背景]:} MCMC方法诞生于1953年Los Alamos国家实验室,Metropolis等人用早期计算机MANIAC研究液态物质状态方程时发明。这是统计物理与计算数学的完美结合,将``用时间换空间''的蒙特卡洛思想与马氏链的收敛理论结合,开创了现代计算统计学。1970年Hastings推广后,成为贝叶斯统计、机器学习、量子场论等领域不可或缺的工具。}

\subsection{马氏链和 MCMC 介绍}

{\color{red}\textbf{[核心挑战]:} 高维积分的``维数灾难''——要计算期望$\mathbb{E}_{\pi}[f(X)] = \int f(x)\pi(x)dx$,如果用网格法,10维空间需要$10^{10}$个点,100维则是$10^{100}$个点(超过宇宙原子数)!传统数值积分在高维空间彻底失效。MCMC的天才之处:不穷举空间,而是设计智能随机游走,自动找到并停留在重要区域。}

实际工作中经常遇到分布复杂的高维随机向量抽样问题。§3.2.4的重要抽样法可以应付维数不太高的情况, 但是对于维数很高而且分布很复杂 (比如, 分布密度多峰而且位置不易确定的情况)则难以处理。

{\color{red}\textbf{[为什么高维困难]:} 高维空间的``维数灾难''——10维空间中要达到1维情况下的采样密度,需要的样本数增长$10^{10}$倍!多峰分布更糟:峰的位置像``大海捞针'',传统方法几乎无法找到所有重要区域。}

MCMC(马氏链蒙特卡洛) 是一种对高维随机向量抽样的方法, 此方法模拟一个马氏链, 使马氏链的平稳分布为目标分布, 由此产生大量的近似服从目标分布的样本, 但样本不是相互独立的。MCMC 的目标分布密度函数或概率函数可以只计算到差一个常数倍的值。MCMC 方法适用范围广, 近年来获得了广泛的应用。

{\color{red}\textbf{[MCMC核心思想]:} 与其试图``一次抽对'',不如构造一个``智能游走''——像醉汉在山地行走,最终会按山的高度分布停留(高处停留概率大)。马氏链保证``无记忆游走'',平稳分布保证``长期收敛到目标''。}

{\color{red}\textbf{[用时间换空间]:} MCMC的哲学是``用时间换空间''——不需要同时存储整个分布(需要指数级内存),只需一个当前状态和转移规则。通过长时间运行马氏链,访问各状态的频率自然收敛到目标分布。这是计算复杂性的根本转变:从空间复杂度指数爆炸,转为时间复杂度多项式增长。}

{\color{red}\textbf{[为什么允许相关样本]:} MCMC产生的样本不独立,但这是可接受的代价!根据遍历定理,相关样本的时间平均仍收敛到真实期望值,只是收敛速度变慢(需乘以``有效样本量因子'')。关键是:相关样本可计算,独立样本不可得——宁要能算的相关,不要算不出的独立。}

先介绍马氏链的概念。

设 \(\left\{  {{X}_{t},t = 0,1,\ldots }\right\}\) 为随机变量序列,称为一个随机过程。称 \({X}_{t}\) 为 “系统在时刻 \(t\) 的状态”。为讨论简单起见,设所有 \({X}_{t}\) 均取值于有限集合 \(S = \{ 1,2,\ldots ,m\}\) ,称 \(S\) 为状态空间。 如果 \(\left\{  {X}_{t}\right\}\) 满足

\[
P\left( {{X}_{t + 1} = j \mid  {X}_{0} = {k}_{0},\ldots {X}_{t - 1} = {k}_{t - 1},{X}_{t} = i}\right)
\]

\[
= P\left( {{X}_{t + 1} = j \mid  {X}_{t} = i}\right)  = {p}_{ij},t = 0,1,\ldots ,{k}_{0},\ldots ,{k}_{t - 1},i,j \in  S, \tag{3.85}
\]

{\color{red}\textbf{[马氏性的物理图像]:} ``无记忆性''——系统的未来只取决于现在,与过去如何到达现在无关。就像抛硬币,下一次正面的概率不管你前面抛了什么。这是最简单的随机演化,却足以描述复杂的平衡分布!}

则称 \(\left\{  {X}_{t}\right\}\) 为马氏链, \({p}_{ij}\) 为转移概率,矩阵 \(P = {\left( {p}_{ij}\right) }_{m \times  m}\) 为转移概率矩阵。显然 \(\mathop{\sum }\limits_{{j = 1}}^{m}{p}_{ij} =\)  \(1,i = 1,2,\ldots ,m\) 。

{\color{red}\textbf{[转移矩阵的物理意义]:} 转移概率矩阵$P$是马氏链的``DNA''——完全编码了系统的动力学规则。$p_{ij}$表示从状态$i$一步跳到状态$j$的概率。行和为1保证概率守恒:从任何状态出发,必然转移到某个状态(包括自身)。矩阵乘法$P^k$给出$k$步转移概率,这是Chapman-Kolmogorov方程的矩阵形式。}

对马氏链, \(P\left( {{X}_{t + k} = j \mid  {X}_{t} = i}\right) \overset{\bigtriangleup }{ = }{p}_{ij}^{\left( k\right) }\) 也不依赖于 \(t\) ,称为 \(k\) 步转移概率。 如果对任意 \(i,j \in  S,i \neq  j\) 都存在 \(k \geq  1\) 使得 \({p}_{ij}^{\left( k\right) } > 0\) 则称 \(\left\{  {X}_{t}\right\}\) 为不可约马氏链。

{\color{red}\textbf{[不可约性的几何图像]:} 不可约性意味着状态空间``强连通''——从任何状态出发,经有限步可达任何其他状态。想象城市地铁网:不可约意味着任意两站都有换乘路径连接。这保证了马氏链能探索整个状态空间,不会永久困在某个子集。对MCMC至关重要:否则链会遗漏部分目标分布!}

不可约马氏链的所有状态是互相连通的,即总能经过若干步后互相转移。对马氏链 \(\left\{  {X}_{t}\right\}\) 的某个状态 \(i\) ,如果存在 \(k \geq  0\) 使得 \({p}_{ii}^{\left( k\right) } > 0\) 并且 \({p}_{ii}^{\left( k + 1\right) } > 0\) ,则称 \(i\) 是非周期的。

{\color{red}\textbf{[非周期性的反例]:} 周期性的直观例子:在棋盘上只允许对角移动的随机游走。从黑格出发,奇数步必在白格,偶数步必在黑格——永远无法在任意时刻达到任意格子。非周期性破坏这种``定时炸弹''模式,允许在任意足够大的时刻返回任何状态,这对收敛到平稳分布是必要的。}

如果一个马氏链所有状态都是非周期的, 则该马氏链称为非周期的。不可约马氏链只要有一个状态是非周期的则所有状态是非周期的。对只有有限个状态的非周期不可约马氏链有

\[
\mathop{\lim }\limits_{{n \rightarrow  \infty }}P\left( {{X}_{n} = j}\right)  = {\pi }_{j},j = 1,2,\ldots ,m, \tag{3.86}
\]

{\color{red}\textbf{[极限定理的深刻性]:} 这是马氏链理论的核心定理!无论初始分布如何,只要满足不可约+非周期,长期行为必然收敛到唯一的极限分布$\pi$。这是``遍历性''的数学表达:系统会``遗忘''初始条件。类比物理:无论初始温度分布如何,孤立系统最终都达到热平衡态——这里$\pi$就是``统计平衡态''。}

其中 \(\left\{  {{\pi }_{j},j = 1,2,\ldots ,m}\right\}\) 为常数,称为 \(\left\{  {X}_{t}\right\}\) 的极限分布。 \(\left\{  {\pi }_{j}\right\}\) 满足方程组

\[
\left\{  \begin{array}{l} \mathop{\sum }\limits_{{i = 1}}^{m}{\pi }_{i}{p}_{ij} = {\pi }_{j},j = 1,2,\ldots ,m \\  \mathop{\sum }\limits_{{j = 1}}^{m}{\pi }_{j} = 1, \end{array}\right.  \tag{3.87}
\]

{\color{red}\textbf{[平稳分布方程的含义]:} 第一个方程$\sum_i \pi_i p_{ij} = \pi_j$是``流入=流出''的平衡条件。用矩阵语言:$\pi P = \pi$,即$\pi$是转移矩阵$P$关于特征值1的左特征向量!这是Perron-Frobenius定理的体现。物理解释:流入状态$j$的总概率流(从所有$i$加权求和)等于$j$本身的概率——动态平衡。第二个方程是归一化条件。}

称满足(3.87)的分布 \(\left\{  {\pi }_{j}\right\}\) 为平稳分布。对只有有限个状态的非周期不可约马氏链,极限分布和平稳分布存在且为同一分布。

{\color{red}\textbf{[为什么叫``平稳'']:} 如果初始分布恰好是$\pi$,即$P(X_0=j)=\pi_j$,那么所有未来时刻的分布都保持为$\pi$:$P(X_t=j)=\pi_j,\forall t$。这就像物理系统的``稳态''——处于平衡态的系统,宏观性质不随时间变化(尽管微观仍在涨落)。对MCMC:平稳分布就是我们想采样的目标分布!}

{\color{red}\textbf{[极限分布=平稳分布的重要性]:} 这个等价性保证了MCMC的可行性!意味着:(1)长时间运行后,链的分布自动收敛到平稳分布,与初值无关;(2)一旦收敛,分布保持不变;(3)时间平均等于系综平均(遍历定理)。这三条是MCMC有效性的数学基础。}

如果允许状态空间 \(S\) 为可列个元素,极限分布的条件需要更多的讨论。对状态 \(i\) ,如从状态 \(i\) 出发总能再返回状态 \(i\) ,则称状态 \(i\) 是常返的 (recurrent)。对常返状态 \(i\) ,如果从 \(i\) 出发首次返回 \(i\) 的时间的期望有限,称 \(i\) 是正常返的。非周期正常返状态称为遍历的。所有状态都遍历的马氏链称为遍历马氏链。非周期遍历马氏链存在唯一的极限分布和平稳分布, 且二者相同。

如果存在 \(\left\{  {{\pi }_{j},j = 1,2,\ldots ,m}\right\}  ,{\pi }_{j} \geq  0,\mathop{\sum }\limits_{{j = 1}}^{m}{\pi }_{j} = 1\) ,使得

\[
{\pi }_{i}{p}_{ij} = {\pi }_{j}{p}_{ji},\forall i \neq  j, \tag{3.88}
\]

{\color{red}\textbf{[细致平衡的深刻含义]:} 这是``微观可逆性''——从$i$到$j$的概率流$\pi_i p_{ij}$与从$j$到$i$的概率流完全相等,系统达到动态平衡。类比:水池进水量=出水量。这比全局平衡(3.87)更强,但更容易构造和验证!}

{\color{red}\textbf{[为什么细致平衡有用]:} 细致平衡是构造MCMC算法的关键工具!它是平稳分布的充分条件(但非必要):满足细致平衡必然满足全局平衡。证明很简单:对(3.88)关于$i$求和,$\sum_i \pi_i p_{ij} = \sum_i \pi_j p_{ji} = \pi_j \sum_i p_{ji} = \pi_j$(利用$\sum_i p_{ji}=1$)。反方向不成立:有净环流的不可逆过程满足全局平衡但不满足细致平衡。}

{\color{red}\textbf{[物理类比——热力学平衡]:} 细致平衡对应物理中的``热力学平衡态''。在平衡态,前向反应率=后向反应率(如$A\rightleftharpoons B$)。这是时间反演对称性的体现:录像正放或倒放都符合物理定律。不满足细致平衡的系统是``耗散系统''(如化学振荡反应),有净能量/概率流动,对应非平衡态。MCMC通常构造满足细致平衡的链,因为更简单。}

称这样的马氏链为细致平衡的 (detailed balanced),这时 \(\left\{  {\pi }_{j}\right\}\) 是 \(\left\{  {X}_{t}\right\}\) 的平稳分布。事实上, 若 \(P\left( {{X}_{t} = i}\right)  = {\pi }_{i},i \in  S\) ,则

\[
P\left( {{X}_{t + 1} = j}\right)  = \mathop{\sum }\limits_{i}{\pi }_{i}{p}_{ij} = \mathop{\sum }\limits_{i}{\pi }_{j}{p}_{ji} = {\pi }_{j}\mathop{\sum }\limits_{i}{p}_{ji} = {\pi }_{j},\forall j \in  S. \tag{3.89}
\]

{\color{red}\textbf{[细致平衡蕴含平稳分布的证明]:} 这个简洁的推导揭示了细致平衡的威力!关键步骤:(1)用细致平衡$\pi_i p_{ij} = \pi_j p_{ji}$,将$\sum_i \pi_i p_{ij}$改写为$\sum_i \pi_j p_{ji}$;(2)提出常数$\pi_j$;(3)利用转移概率归一化$\sum_i p_{ji}=1$。三步证明完毕!这说明设计MCMC时,只需验证细致平衡,自动保证平稳分布。}

马氏链的概念可以推广到 \({X}_{t}\) 的取值集合 \(\mathcal{X}\) 为可列集或 \({\mathbb{R}}^{d}\) 的区域的情形。如果各 \(\left\{  {X}_{t}\right\}\) 的有限维分布是连续型的,则 (3.85) 可以改用条件密度表示,这时的 \(\left\{  {X}_{t}\right\}\) 按照随机过程论中的习惯应该称作马氏过程, 但这里还是叫做马氏链。

如果遍历的不可约马氏链 \(\left\{  {X}_{t}\right\}\) 有平稳分布 \(\pi \left( x\right) ,x \in  \mathcal{X}\) ,则从任意初值出发模拟产生序列 \(\left\{  {X}_{t}\right\}\) ,当 \(t\) 很大时, \({X}_{t}\) 的分布就近似服从 \(\pi\) ,抛弃开始的一段后的 \({X}_{t}\) 序列可以作为分布 \(\pi\) 的相关的样本,抛弃的一段序列叫做老化期。

{\color{red}\textbf{[老化期(Burn-in)的必要性]:} 为什么需要丢弃初始样本?因为初始分布通常远离目标分布$\pi$。想象一个在山谷底部开始的随机游走,目标分布集中在山顶——初期样本严重偏向谷底,不能代表$\pi$!老化期让链``爬到''高概率区域。但多长算够?没有通用公式,需要诊断:观察统计量(如均值)的轨迹,当其稳定在某个值附近震荡,说明接近收敛。}

{\color{red}\textbf{[如何选择初值]:} 初值选择策略:(1)随机选取——最保守,但可能需要很长老化期;(2)利用先验知识——在高概率区域附近初始化(如用最大似然估计或MAP估计);(3)过度分散初始化(Overdispersed Initialization)——多条链从完全不同的极端位置开始,若都收敛到同一分布,增强信心。好的初值能显著减少老化期,但不影响最终结果(仅影响效率)。}

设 \(Y \sim  \pi \left( \cdot \right) ,h\left( y\right) ,y \in  \mathcal{X}\) 是有界函数,为估计 \(\theta  \triangleq  {Eh}\left( Y\right)\) ,用 \({X}_{t},t = k + 1,\ldots ,n\) 作为 \(\pi\) 的样本,用估计量 \(\widehat{\theta } = \frac{1}{n - k}\mathop{\sum }\limits_{{t = k + 1}}^{n}h\left( {X}_{t}\right)\) 来估计 \(\theta\) ,则 \(\widehat{\theta }\) 是 \(\theta\) 的强相合估计。老化期长度 \(k\) 可以从 \(\widehat{\theta }\) 的变化图形经验地选取。

{\color{red}\textbf{[遍历定理保证收敛]:} $\widehat{\theta}$的强相合性来自遍历定理:对遍历马氏链,时间平均几乎必然收敛到空间平均,即$\frac{1}{n}\sum_{t=1}^n h(X_t) \xrightarrow{a.s.} \mathbb{E}_\pi[h(X)]$。这是大数定律在相关样本上的推广!虽然样本相关,但由于链``遍历''整个空间,长期平均仍然正确。这是MCMC方法论的核心数学保证。}

这样的估计量 \(\widehat{\theta }\) 是相关样本的平均值,无法用原来独立样本的公式估计 \(\operatorname{Var}\left( \widehat{\theta }\right)\) 从而得到 \(\widehat{\theta }\) 的标准误差。

{\color{red}\textbf{[有效样本量]:} 相关样本的方差估计需要考虑自相关!对独立样本,$\text{Var}(\bar{X}) = \sigma^2/n$。但对相关样本,真实方差是$\text{Var}(\widehat{\theta}) = \frac{\sigma^2}{n}(1 + 2\sum_{k=1}^\infty \rho_k)$,其中$\rho_k$是滞后$k$的自相关系数。定义有效样本量$n_{\text{eff}} = n/(1+2\sum \rho_k)$——这是``等价独立样本数''。如果自相关很强,$n_{\text{eff}}\ll n$,说明样本信息量严重打折!}

为了估计 \(\operatorname{Var}\left( \widehat{\theta }\right)\) ,可以采用如下的分段平均法。把样本 \({X}_{k + 1},\ldots ,{X}_{n}\) 分为 \(s\) 段,每段 \(r\) 个 (设 \(n - k = {sr}\) )。设第 \(j\) 段的 \(r\) 个 \(h\left( {X}_{t}\right)\) 的平均值为 \({Z}_{j},j = 1,2,\ldots ,s\) ,设 \(\left\{  {{Z}_{j},j = 1,2,\ldots ,s}\right\}\) 的样本方差为 \({\widehat{\sigma }}^{2} = \frac{1}{s - 1}\mathop{\sum }\limits_{{j = 1}}^{s}{\left( {Z}_{j} - \bar{Z}\right) }^{2}\) ,因为 \(\widehat{\theta }\) 等于 \(\left\{  {{Z}_{j},j = 1,2,\ldots ,s}\right\}\) 的样本均值 \(\bar{Z}\) ,当 \(r\) 足够大时,可以认为各 \({Z}_{j}\) 相关性已经很弱,这时 \(\widehat{\theta }\) 的方差可以用 \({\widehat{\sigma }}^{2}/s\) 估计。 \(r\) 的大小依赖于不同时刻的 \({X}_{t}\) 的相关性强弱,相关性越强,需要的 \(r\) 越大。

{\color{red}\textbf{[批处理方法的直觉]:} 分段平均法(Batch Means)的思想:将长链分成$s$段,每段$r$个样本。若$r$足够大(大于自相关长度),段内样本的平均值$Z_j$之间近似独立!这样就把相关样本问题转化为独立样本问题。类比:测量强相关的温度序列,如果每隔足够长时间取一次平均,这些平均值之间就近似独立了。选择$r$的经验法则:观察自相关图,当$\rho_r < 0.05$时的$r$。}

以上的方法就是 MCMC 方法 (马氏链蒙特卡洛)。一般地,对高维或取值空间 \(\mathcal{X}\) 结构复杂的随机向量 \(X\) , MCMC 方法构造取值于 \(\mathcal{X}\) 的马氏链,使其平稳分布为 \(X\) 的目标分布。 模拟此马氏链,抛弃开始的部分抽样值,把剩余部分作为 \(X\) 的非独立抽样。非独立抽样的估计效率比独立抽样低。

MCMC 方法的关键在于如何从第 \(t\) 时刻转移到第 \(t + 1\) 时刻。好的转移算法应该使得马氏链比较快地收敛到平稳分布,并且不会长时间地停留在取值空间 \(\mathcal{X}\) 的局部区域内(在目标分布是多峰分布且峰高度差异较大时容易出现这种问题)。

{\color{red}\textbf{[好的MCMC算法的标准]:} 评价MCMC算法的两大指标:(1)\textbf{收敛速度}——到达平稳分布需要多少步?越快越好,减少老化期浪费。(2)\textbf{混合性质(Mixing)}——能否高效探索整个目标分布?差的混合导致链``卡住''在局部区域(如低概率峰),遗漏重要区域。理想算法:快速收敛+良好混合。这两者往往矛盾:大步跳跃利于混合但降低接受率,小步提高接受率但减慢探索。}

{\color{red}\textbf{[多峰分布的挑战]:} 多峰分布是MCMC的``噩梦''!想象两座高山被深谷分隔:如果链困在一座山上,很难跨越低概率谷底到达另一座山。标准MCMC倾向于停留在最高峰附近,遗漏次高峰。后果:估计严重偏差!解决方案:(1)平行回火(Parallel Tempering)——同时运行多条不同``温度''的链,高温链更容易跨越峡谷;(2)自适应步长——动态调整提议分布;(3)Hamiltonian Monte Carlo——利用梯度信息导航。}

Metropolis-Hasting 方法 (MH 方法) 是一个基本的 MCMC 算法, 此算法在每一步试探地进行转移 (如随机游动),如果转移后能提高状态 \({x}_{t}\) 在目标分布 \(\pi\) 中的密度值则接受转移结果, 否则以一定的概率决定是转移还是停留不动。

{\color{red}\textbf{[历史背景]:} Metropolis等人1953年在Los Alamos用MANIAC计算机研究液态物质方程,发明了这个算法。Hastings 1970年推广。这是统计物理与计算统计的完美结合——用``蒙特卡洛''模拟平衡态分布!}

{\color{red}\textbf{[Metropolis-Hastings的天才设计]:} MH算法的核心是``提议-修正''两步走:(1)\textbf{提议}——从当前点$x$,用简单的提议分布$T(y|x)$(如正态分布)随机生成候选点$y$。提议可以很``愚蠢'',甚至与目标分布无关!(2)\textbf{修正}——用精巧设计的接受概率$r(x,y)$,修正提议的偏差,强制满足细致平衡。这种分离使算法极其灵活:提议负责探索,接受概率负责正确性。}

Gibbs 抽样是另外一种常用的 MCMC 方法, 此方法轮流延各个坐标轴方向转移, 且转移概率由当前状态下用其它坐标预测转移方向坐标的条件分布给出。因为利用了目标分布的条件分布,所以 Gibbs 抽样方法的效率比 MH 方法效率更高。

\subsection{Metropolis-Hasting 抽样}

设随机变量 \(X\) 分布为 \(\pi \left( x\right) ,x \in  \mathcal{X}\) 。为论述简单起见仍假设 \(\mathcal{X}\) 是离散集合。算法需要一个试转移概率函数 \(T\left( {y \mid  x}\right) ,x,y \in  \mathcal{X}\) ,满足 \(0 \leq  T\left( {y \mid  x}\right)  \leq  1,\mathop{\sum }\limits_{y}T\left( {y \mid  x}\right)  = 1\) ,并且

\[
T\left( {y \mid  x}\right)  > 0 \Leftrightarrow  T\left( {x \mid  y}\right)  > 0. \tag{3.90}
\]

算法首先从 \(\mathcal{X}\) 中任意取初值 \({X}^{\left( 0\right) }\) 。设经过 \(t\) 步后算法的当前状态为 \({X}^{\left( t\right) }\) ,则下一步由试转移分布 \(T\left( {y \mid  {X}^{\left( t\right) }}\right)\) 抽取 \(Y\) ,并生成 \(U \sim  \mathrm{U}\left( {0,1}\right)\) ,然后按如下规则转移:

\[
{X}^{\left( t + 1\right) } = \left\{  \begin{array}{ll} Y & \text{ 若 }U \leq  r\left( {{X}^{\left( t\right) },Y}\right) \\  {X}^{\left( t\right) } & \text{ 否则 } \end{array}\right.  \tag{3.91}
\]

其中

\[
r\left( {x,y}\right)  = \min \left\{  {1,\frac{\pi \left( y\right) T\left( {x \mid  y}\right) }{\pi \left( x\right) T\left( {y \mid  x}\right) }}\right\}  . \tag{3.92}
\]

{\color{red}\textbf{[接受率的精妙设计]:} 为什么这样设计就能工作?关键是``向上爬山总接受,向下走按比例接受''。如果新点$y$比当前点$x$的$\pi$值大,必定接受(爬高山);如果更小,按比例$\pi(y)/\pi(x)$接受(偶尔下坡避免卡住)。这保证了细致平衡!}

{\color{red}\textbf{[接受率公式的两部分解读]:} 接受率$r(x,y) = \min\{1, \frac{\pi(y)T(x|y)}{\pi(x)T(y|x)}\}$包含两个比值:(1)$\frac{\pi(y)}{\pi(x)}$——目标分布的比值,衡量新点$y$相对于当前点$x$的``重要性'';(2)$\frac{T(x|y)}{T(y|x)}$——提议分布的比值,修正提议的不对称性。如果$T$对称(如$T(y|x)=T(x|y)$),第二项为1,简化为Metropolis算法!这个公式保证:即使提议分布``偏心''(如更容易向某方向跳),接受率会自动纠偏。}

{\color{red}\textbf{[接受率的物理直觉——Boltzmann分布]:} 在统计物理中,若$\pi(x) \propto e^{-E(x)/T}$(Boltzmann分布,$E(x)$是能量,$T$是温度),则$\frac{\pi(y)}{\pi(x)} = e^{-(E(y)-E(x))/T} = e^{-\Delta E/T}$。接受率变为$\min\{1, e^{-\Delta E/T}\}$。解释:(1)若$\Delta E < 0$(能量降低),必然接受——系统自发趋向低能态;(2)若$\Delta E > 0$(能量升高),以概率$e^{-\Delta E/T}$接受——这是热涨落!温度$T$越高,越容易接受``不利''移动,帮助逃离局部极小值。这就是模拟退火的原理。}

在 MH 算法中如果取 \(T\left( {y \mid  x}\right)  = T\left( {x \mid  y}\right)\) ,则 \(r\left( {x,y}\right)  = \min \left( {1,\frac{\pi \left( y\right) }{\pi \left( x\right) }}\right)\) ,相应的算法称为 Metropolis 抽样法。

{\color{red}\textbf{[Metropolis算法的简化]:} 对称提议分布的常见例子:$T(y|x) = \mathcal{N}(y; x, \sigma^2 I)$——以当前点为中心的各向同性正态分布。从$x$提议到$y$的概率等于从$y$提议到$x$的概率,故对称。此时接受率极度简化为$\min\{1, \pi(y)/\pi(x)\}$,只需计算目标分布的比值!这是最常用的MCMC变体,因其实现简单且理论清晰。}

如果取 \(T\left( {y \mid  x}\right)  = g\left( y\right)\) (不依赖于 \(x\) ),则 \(r\left( {x,y}\right)  = \min \left( {1,\frac{\pi \left( y\right) /g\left( y\right) }{\pi \left( x\right) /g\left( x\right) }}\right)\) ,相应的算法称为 Metropolis 独立抽样法,和重要抽样有相似之处,试抽样分布 \(g\left( \cdot \right)\) 经常取为相对重尾的分布。

{\color{red}\textbf{[独立抽样的特殊性]:} 独立MH的提议$T(y|x) = g(y)$与当前状态$x$无关——每次从同一分布$g$独立抽样!这看似违反``马氏链''精神,但接受/拒绝步骤引入了对$x$的依赖。接受率$\min\{1, \frac{\pi(y)/g(y)}{\pi(x)/g(x)}\}$中的$\pi/g$是重要性权重!这与重要性采样的联系:如果$g$接近$\pi$,权重比接近1,高接受率;如果$g$偏离$\pi$,低接受率。选择$g$的关键:重尾(覆盖$\pi$的支撑),但不能太平坦(导致低效拒绝)。}

在 MH 算法中,目标分布 \(\pi \left( x\right)\) 可以用差一个常数倍的 \(\widetilde{\pi }\left( x\right)  = {C\pi }\left( x\right)\) 代替,这样关于目标分布仅知道差一个常数倍的 \(\widetilde{\pi }\left( x\right)\) 的情形,也可以使用此算法。

{\color{red}\textbf{[为什么不需要归一化常数]:} 因为接受率$r(x,y)$中的常数$C$会上下抵消!这是MCMC的巨大优势——贝叶斯后验$p(\theta|{\rm data})\propto p({\rm data}|\theta)p(\theta)$不需要计算积分$\int p({\rm data}|\theta)p(\theta)d\theta$,直接用未归一化的即可。}

{\color{red}\textbf{[归一化常数消除的数学]:} 若$\pi(x) = \widetilde{\pi}(x)/Z$,其中$Z = \int \widetilde{\pi}(x)dx$未知,则$\frac{\pi(y)}{\pi(x)} = \frac{\widetilde{\pi}(y)/Z}{\widetilde{\pi}(x)/Z} = \frac{\widetilde{\pi}(y)}{\widetilde{\pi}(x)}$,$Z$完美消去!这对贝叶斯推断是革命性的:后验$p(\theta|D) = \frac{p(D|\theta)p(\theta)}{p(D)}$,其中证据$p(D) = \int p(D|\theta)p(\theta)d\theta$通常不可计算(高维积分!)。MH只需$\widetilde{\pi}(\theta) = p(D|\theta)p(\theta)$(似然×先验),无需归一化!这打开了复杂贝叶斯模型的大门。}

下面说明 \(\mathrm{{MH}}\) 抽样方法的合理性。我们来验证 \(\mathrm{{MH}}\) 抽样的转移概率 \(A\left( {x,y}\right)  = P\left( {{X}^{\left( t + 1\right) } = }\right.\)  \(\left. {y \mid  {X}^{\left( t\right) } = x}\right)\) 满足细致平衡条件。

{\color{red}\textbf{[为什么需要证明细致平衡]:} MH算法的正确性完全依赖于其满足细致平衡条件$\pi(x)A(x,y) = \pi(y)A(y,x)$。这不是显然的!因为算法包含两个随机步骤(提议+接受/拒绝),它们的组合效果必须精确满足这个平衡关系。下面的证明揭示接受率$r(x,y)$的精巧设计如何保证此性质——这是MH算法天才之处的数学基础。}

易见

\[
A\left( {x,y}\right)  = \left\{  \begin{array}{ll} T\left( {y \mid  x}\right) r\left( {x,y}\right) , & y \neq  x, \\  T\left( {x \mid  x}\right)  + \mathop{\sum }\limits_{{z \neq  x}}T\left( {z \mid  x}\right) \left\lbrack  {1 - r\left( {x,z}\right) }\right\rbrack  , & y = x, \end{array}\right.  \tag{3.93}
\]

{\color{red}\textbf{[转移概率的两部分]:} MH的转移概率$A(x,y)$有两种情况:(1)$y \neq x$:必须先提议到$y$(概率$T(y|x)$),再被接受(概率$r(x,y)$),故$A(x,y) = T(y|x)r(x,y)$;(2)$y=x$(停留):两种途径——直接提议自己(概率$T(x|x)$,通常为0),或提议其他点$z$但被拒绝(概率$T(z|x)[1-r(x,z)]$)。第二项是所有被拒绝提议的总和,保证$\sum_y A(x,y) = 1$(概率归一化)。}

于是当 \(x \neq  y\) 时

\[
\pi \left( x\right) A\left( {x,y}\right)  = \pi \left( x\right) T\left( {y \mid  x}\right) \min \left\{  {1,\frac{\pi \left( y\right) T\left( {x \mid  y}\right) }{\pi \left( x\right) T\left( {y \mid  x}\right) }}\right\}   = \min \{ \pi \left( x\right) T\left( {y \mid  x}\right) ,\pi \left( y\right) T\left( {x \mid  y}\right) \} ,
\]

{\color{red}\textbf{[证明的关键步骤]:} 这一步是证明的核心技巧!将$\pi(x)$乘进$\min$内部:$\pi(x) \min\{1, \frac{\pi(y)T(x|y)}{\pi(x)T(y|x)}\} = \min\{\pi(x), \frac{\pi(y)T(x|y)}{T(y|x)}\} = \min\{\pi(x)T(y|x), \pi(y)T(x|y)\}$(利用$\min\{a, bc/a\} = \min\{a, bc\}/\max\{1, a/bc\}$的性质)。关键观察:最终表达式对$x$和$y$是\textbf{对称的}!}

等式右侧关于 \(x,y\) 是对称的,所以等式左侧把 \(x,y\) 交换后仍相等。所以, \(\mathrm{{MH}}\) 构造的马氏链以 \(\{ \pi \left( x\right) \}\) 为平稳分布。多数情况下 MH 构造的马氏链也以 \(\{ \pi \left( x\right) \}\) 为极限分布。

{\color{red}\textbf{[证明完成]:} $\min\{\pi(x)T(y|x), \pi(y)T(x|y)\}$关于$(x,y)$对称,意味着$\pi(x)A(x,y) = \pi(y)A(y,x)$——细致平衡成立!结论:无论提议分布$T$是什么,只要按MH公式设计接受率$r(x,y)$,就自动满足细致平衡,从而$\pi$是平稳分布。这个证明展示了MH算法的优美性:用一个简单的代数技巧(对称化),将任意提议分布转化为正确的采样器。}

(3.92) 中的 \(r\left( {x,y}\right)\) 还可以推广为如下的形式

\[
\widetilde{r}\left( {x,y}\right)  = \frac{\alpha \left( {x,y}\right) }{\pi \left( x\right) T\left( {y \mid  x}\right) }, \tag{3.94}
\]

其中 \(\alpha \left( {x,y}\right)\) 是任意的满足 \(\alpha \left( {x,y}\right)  = \alpha \left( {y,x}\right)\) 且使得 \(\widetilde{r}\left( {x,y}\right)  \leq  1\) 的函数。易见这样的 \(\widetilde{r}\left( {x,y}\right)\) 仍使得生成的马氏链满足细致平衡条件。

例 3.7.1. \(X\) 的取值集合 \(\mathcal{X}\) 可能是很大的,以至于无法穷举,目标分布 \(\pi \left( x\right)\) 可能是只能确定到差一个常数倍。

例如, 设

\[
\mathcal{X} = \left\{  {\mathbf{x} = \left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)  : \left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right) \text{ 为 }\left( {1,2,\ldots ,n}\right) \text{ 的一个排列,并满足 }\mathop{\sum }\limits_{{j = 1}}^{n}j{x}_{j} > a}\right\}  ,
\]

其中 \(a\) 是一个给定的常数。用 \(\left| \mathcal{X}\right|\) 表示 \(\mathcal{X}\) 的元素个数,当 \(n\) 较大时 \(\mathcal{X}\) 是 \(\left( {1,2,\ldots ,n}\right)\) 的所有 \(n!\) 个排列的一个子集, \(\left| \mathcal{X}\right|\) 很大,很难穷举 \(\mathcal{X}\) 的元素,从而 \(\left| \mathcal{X}\right|\) 未知。

设 \(X\) 服从 \(\mathcal{X}\) 上的均匀分布,即 \(\pi \left( \mathbf{x}\right)  = C,x \in  \mathcal{X},C = 1/\left| \mathcal{X}\right|\) 但 \(C\) 未知。要用 MH 方法产生 \(X\) 的抽样序列。

试抽样 \(T\left( {\mathbf{y} \mid  \mathbf{x}}\right)\) 如果允许转移到所有的 \(\mathbf{y}\) 是很难执行的,因为 \(\mathbf{y}\) 的个数太多了。我们定义一个 \(\mathbf{x}\) 的近邻的概念,仅考虑转移到 \(\mathbf{x}\) 的近邻。一种定义是,如果把 \(\mathbf{x}\) 的 \(n\) 个元素中的某两个交换位置后可以得到 \(\mathbf{y} \in  \mathcal{X}\) ,则 \(\mathbf{y}\) 称为 \(\mathbf{x}\) 的一个近邻,记 \(N\left( \mathbf{x}\right)\) 为 \(\mathbf{x}\) 的所有近邻的集合,记 \(\left| {N\left( \mathbf{x}\right) }\right|\) 为 \(\mathbf{x}\) 的近邻的个数。当 \(n\) 很大时,求 \(N\left( \mathbf{x}\right)\) 也需要从 \({C}_{n}^{2} = \frac{1}{2}n\left( {n - 1}\right)\) 个可能的元素中用穷举法选择。取试转移概率函数为

\[
T\left( {\mathbf{y} \mid  \mathbf{x}}\right)  = \frac{1}{\left| N\left( \mathbf{x}\right) \right| },\mathbf{x},\mathbf{y} \in  \mathcal{X},
\]

即从 \(x\) 出发,等可能地试转移到 \(x\) 的任何一个近邻上。

因为目标分布 \(\pi \left( \mathbf{x}\right)\) 是常数,所以这时

\[
r\left( {\mathbf{x},\mathbf{y}}\right)  = \min \left( {1,\frac{\left| N\left( \mathbf{x}\right) \right| }{\left| N\left( \mathbf{y}\right) \right| }}\right) ,
\]

即从 \(\mathbf{x}\) 试转移到 \(\mathbf{y}\) 后,如果 \(\mathbf{y}\) 的近邻数不超过 \(\mathbf{x}\) 的近邻数则确定转移到 \(\mathbf{y}\) ,否则,仅按概率 \(\left| {N\left( \mathbf{x}\right) }\right| /\left| {N\left( \mathbf{y}\right) }\right|\) 转移到 \(\mathbf{y}\) 。这就构成了对 \(X\) 抽样的 \(\mathrm{{MH}}\) 算法。 连续型分布的 \(\mathbf{{MH}}\) 抽样法 对于连续型的目标分布,设 \(\pi \left( x\right)\) 为目标分布的密度,这时 \(T\left( {y \mid  x}\right)\) 改为给定 \(x\) 条件下的试抽样密度, \(r\left( {x,y}\right)\) 定义不变,算法和离散型目标分布的情形相同。

例 3.7.2. 考虑一个贝叶斯推断问题。在金融投资中, 投资者经常把若干种证券组合在一起来减少风险。假设有 5 支股票的 \(n = {250}\) 个交易日的收益率记录,每个交易日都找出这 5 支股票收益率最高的一个,设 \({X}_{i}\) 表示第 \(i\) 支股票在 \(n\) 个交易日中收益率为最高的次数 \(\left( {i = 1,2,\ldots ,5}\right)\) 。设 \(\left( {{X}_{1},\ldots ,{X}_{5}}\right)\) 服从多项分布,相应的概率假设为

\[
\mathbf{p} = \left( {\frac{1}{3},\frac{1 - \beta }{3},\frac{1 - {2\beta }}{3},\frac{2\beta }{3},\frac{\beta }{3}}\right) ,
\]

其中 \(\beta  \in  \left( {0,{0.5}}\right)\) 为未知参数。假设 \(\beta\) 有先验分布 \({p}_{0}\left( \beta \right)  \sim  \mathrm{U}\left( {0,{0.5}}\right)\) 。设 \(\left( {{x}_{1},\ldots ,{x}_{5}}\right)\) 为 \(\left( {{X}_{1},\ldots ,{X}_{5}}\right)\) 的观测值,则 \(\beta\) 的后验分布为

\[
f\left( {\beta  \mid  {x}_{1},\ldots ,{x}_{5}}\right)  \propto  p\left( {{x}_{1},\ldots ,{x}_{5} \mid  \beta }\right) {p}_{0}\left( \beta \right)
\]

\[
= \left( \begin{matrix} n \\  {x}_{1},\ldots ,{x}_{5} \end{matrix}\right) {\left( \frac{1}{3}\right) }^{{x}_{1}}{\left( \frac{1 - \beta }{3}\right) }^{{x}_{2}}{\left( \frac{1 - {2\beta }}{3}\right) }^{{x}_{3}}{\left( \frac{2\beta }{3}\right) }^{{x}_{4}}{\left( \frac{\beta }{3}\right) }^{{x}_{5}}\frac{1}{0.5}{I}_{\left( 0,{0.5}\right) }\left( \beta \right)
\]

\[
\propto  {\left( 1 - \beta \right) }^{{x}_{2}}{\left( 1 - 2\beta \right) }^{{x}_{3}}{\beta }^{{x}_{4} + {x}_{5}}{I}_{\left( 0,{0.5}\right) }\left( \beta \right) \overset{ \bigtriangleup  }{ = }\widetilde{\pi }\left( \beta \right) .
\]

为了求 \(\beta\) 后验均值,需要产生服从 \(f\left( {\beta  \mid  {x}_{1},\ldots ,{x}_{5}}\right)\) 的抽样。从 \(\beta\) 的后验分布很难直接抽样,采用 Metropolis 抽样法。设当前 \(\beta\) 的状态为 \({\beta }^{\left( t\right) }\) ,取试抽样分布 \(T\left( {y \mid  {\beta }^{\left( t\right) }}\right)\) 为 \(\mathrm{U}\left( {0,{0.5}}\right)\) , 则 \(T\left( {y \mid  x}\right)  = T\left( {x \mid  y}\right)\) ,

\[
r\left( {{\beta }^{\left( t\right) },y}\right)  = \min \left( {1,\frac{\widetilde{\pi }\left( y\right) }{\widetilde{\pi }\left( {\beta }^{\left( t\right) }\right) }}\right)  = \min \left( {1,{\left( \frac{1 - y}{1 - {\beta }^{\left( t\right) }}\right) }^{{x}_{2}}{\left( \frac{1 - {2y}}{1 - 2{\beta }^{\left( t\right) }}\right) }^{{x}_{3}}{\left( \frac{y}{{\beta }^{\left( t\right) }}\right) }^{{x}_{4} + {x}_{5}}}\right) ,
\]

从 \(\mathrm{U}\left( {0,{0.5}}\right)\) 试抽取 \(y\) ,以概率 \(r\left( {{\beta }^{\left( t\right) },y}\right)\) 接受 \({\beta }^{\left( t + 1\right) } = y\) 即可。

随机游动 \(\mathbf{{MH}}\) 算法 \({MH}\) 抽样中试转移概率函数 \(T\left( {y \mid  x}\right)\) 较难找到,容易想到的是从 \({x}^{\left( t\right) }\) 作随机游动的试转移方法, 叫做随机游动 Metropolis 抽样。

{\color{red}\textbf{[类比]:} 随机游动MH就像``醉汉爬山''——每步随机选个方向走一小步,如果上坡就走,如果下坡就``掷骰子''决定。步长$\sigma$是关键:太大则总被拒绝(想跳太远摔下来),太小则蜗牛爬行(几万步还在山脚)。}

{\color{red}\textbf{[随机游走的局部性]:} 随机游走MH是\textbf{局部探索}策略——每步只在当前点附近小范围移动,逐步探索空间。优点:简单,不需要对目标分布的先验知识;缺点:探索慢,在高维空间效率低(需指数级步数才能遍历空间)。对比独立MH的全局跳跃:随机游走更保守但更稳定。这是探索(Exploration)与利用(Exploitation)的权衡。}

设 \(X\) 的目标分布 \(\pi \left( \mathbf{x}\right)\) 取值于欧式空间 \(\mathcal{X} = {\mathbb{R}}^{d}\) 。从 \({\mathbf{x}}^{\left( t\right) }\) 出发试转移,令

\[
\mathbf{y} = {\mathbf{x}}^{\left( t\right) } + {\mathbf{\varepsilon }}_{t}, \tag{3.95}
\]

{\color{red}\textbf{[提议机制]:} 提议$\mathbf{y} = \mathbf{x}^{(t)} + \varepsilon_t$是当前位置加随机扰动——这是Brownian运动的离散化!$\varepsilon_t$的分布$g$决定探索特性。对称性$g(\varepsilon) = g(-\varepsilon)$保证$T(y|x) = g(y-x) = g(x-y) = T(x|y)$,简化接受率为$\min\{1, \pi(y)/\pi(x)\}$。这个对称性至关重要:破坏它需要Hastings修正,增加计算负担。}

其中 \({\varepsilon }_{t} \sim  g\left( {\mathbf{x};\sigma }\right)\) 对不同 \(t\) 是独立同分布的, \(T\left( {\mathbf{y} \mid  \mathbf{x}}\right)  = g\left( {\mathbf{y} - \mathbf{x}}\right)\) 。设 \(g\) 是关于 \(\mathbf{x} = \mathbf{0}\) 对称的分布,则 \(T\left( {\mathbf{y} \mid  \mathbf{x}}\right)  = T\left( {\mathbf{x} \mid  \mathbf{y}}\right)\) 。

常取 \(g\) 为 \(\mathrm{N}\left( {\mathbf{0},{\sigma }^{2}I}\right)\) 和半径为 \(\sigma\) 的中心为0的球内的均匀分布。

{\color{red}\textbf{[两种常见提议分布]:} (1)\textbf{各向同性正态}$\mathcal{N}(0, \sigma^2 I)$——最流行!各方向独立,方差相同。$\sigma$控制步长:$\sigma$小→短步(高接受率,慢探索),$\sigma$大→大步(低接受率,快探索)。(2)\textbf{均匀球}$\text{Unif}(B_\sigma(0))$——在半径$\sigma$球内均匀分布。与正态类似,但有硬边界(正态是软边界)。选择原则:正态更常用,因其中心极限定理保证的良好性质;均匀球在需要严格步长上界时有用。}

转移法则为: 从 \({\mathbf{x}}^{\left( t\right) }\) 出发试转移到 \(\mathbf{y}\) 后,若 \(\pi \left( \mathbf{y}\right)  > \pi \left( {\mathbf{x}}^{\left( t\right) }\right)\) 则令 \({\mathbf{x}}^{\left( t + 1\right) } = \mathbf{y}\) ; 否则,独立地抽取 \(U \sim  \mathrm{U}\left( {0,1}\right)\) ,取

\[
{\mathbf{x}}^{\left( t + 1\right) } = \left\{  \begin{array}{ll} \mathbf{y}, & \text{ 当 }U \leq  \pi \left( \mathbf{y}\right) /\pi \left( {\mathbf{x}}^{\left( t\right) }\right) , \\  {\mathbf{x}}^{\left( t\right) }, & \text{ 否则. } \end{array}\right.
\]

{\color{red}\textbf{[接受/拒绝的直觉]:} 这个规则体现了``贪婪+随机''的平衡:(1)若$\pi(y) > \pi(x)$(爬山),总是接受——利用梯度信息;(2)若$\pi(y) < \pi(x)$(下山),以概率$\pi(y)/\pi(x)$接受——这是关键!下坡概率与密度比成正比,保证长期停留频率正确。类比退火:偶尔接受``坏''移动,避免困在局部最优。若总拒绝下山,链会永远停在首个局部极大值!}

随机游动 \(\mathrm{{MH}}\) 算法是一种 Metropolis 抽样方法。随机游动的步幅 \(\sigma\) 是重要参数,步幅过大导致拒绝率大, 步幅过小使得序列的相关性太强, 收敛到平衡态速度太慢。一个建议选法是试验各种选法, 使得试抽样被接受的概率在 0.25 到 0.35 之间。

{\color{red}\textbf{[步长调优的理论]:} 接受率与步长$\sigma$的关系:$\sigma \to 0$时,接受率→1(总是接受),但步长趋零,探索极慢;$\sigma \to \infty$时,接受率→0(总被拒绝),链几乎不动。最优在中间!Roberts等人(1997)证明:对高维正态目标分布,\textbf{最优接受率约为0.234}(23.4\%规则)。这是一个惊人的普适结果!虽然针对特定情况,但对许多复杂分布也适用。实践中,0.25-0.35是经验黄金区间。}

{\color{red}\textbf{[自适应调节策略]:} 实际应用中,常在老化期\textbf{自适应调整}$\sigma$:每隔固定步数,检查最近的接受率;若过高(>0.5),增大$\sigma$(如乘以1.1);若过低(<0.2),减小$\sigma$(如除以1.1)。迭代直到接受率稳定在目标区间。\textbf{注意}:调整必须在老化期完成,收敛后不能再调(会破坏遍历性)!这个技巧极大提高实用性,避免手动试错。}

例 3.7.3. 考虑如下的简单气体模型: 在平面区域 \(G = \left\lbrack  {0,A}\right\rbrack   \times  \left\lbrack  {0,B}\right\rbrack\) 内有 \(K\) 个直径为 \(d\) 的刚性圆盘。随机向量 \(\mathbf{X} = \left( {{x}_{1},{y}_{1},\ldots ,{x}_{K},{y}_{k}}\right)\) 为这些圆盘的位置坐标。分布 \(\pi \left( \mathbf{x}\right)\) 是 \(G\) 内所有允许位置的均匀分布。希望对 \(\pi\) 抽样。

先找一个初始的允许位置 \({\mathbf{x}}^{\left( 0\right) }\) 。比如,把圆盘整齐地排列在左上角。

设已得到 \({\mathbf{x}}^{\left( t\right) }\) ,随机选取一个圆盘 \(i\) ,把圆盘 \(i\) 的位置试移动到 \(\left( {{x}_{i}^{\prime },{y}_{i}^{\prime }}\right)  = \left( {{x}_{i} + {\delta }_{i},{y}_{i} + {\epsilon }_{i}}\right)\) , 其中 \({\delta }_{i},{\epsilon }_{i}\) 独立同 \(\mathrm{N}\left( {0,{\sigma }^{2}}\right)\) 分布。如果得到的位置是允许的则接受结果,否则留在原地不动。

\subsection{Gibbs 抽样}

一般的 MH 抽样每一步首先进行尝试运动, 然后根据新的状态是否靠近目标分布来接受或拒绝试抽样点, 所以可能会存在多次的无效尝试, 效率较低。

{\color{red}\textbf{[MH的拒绝损失]:} MH算法的``接受/拒绝''机制虽然保证正确性,但带来效率损失:每次拒绝意味着浪费一次计算(计算了$\pi(y)$但没采用)。高拒绝率→低有效样本量。在高维空间或步长不当时,拒绝率可能高达90\%!能否设计``总是接受''的MCMC?Gibbs抽样给出肯定答案——代价是需要条件分布可采样。}

Gibbs 抽样是另外一种 MCMC 方法, 它仅在坐标轴方向尝试转移, 用当前点的条件分布决定下一步的试抽样分布, 所有试抽样都被接受, 不需要拒绝, 所以效率可以更高。

{\color{red}\textbf{[Gibbs抽样的巧妙之处]:} 为什么从条件分布抽样就不需要拒绝?因为条件分布$p(x_i|\mathbf{x}_{(-i)})$本身就是目标分布的``切片''!就像修房子,MH是``整体平移试探'',Gibbs是``每次只修一个房间,按正确图纸修'',当然不会出错。}

{\color{red}\textbf{[Gibbs vs MH的对比]:} \textbf{Gibbs优势}:(1)无拒绝,接受率100\%;(2)自动适应目标分布形状(通过条件分布);(3)计算效率高(每步必然移动)。\textbf{Gibbs劣势}:(1)需要知道并能采样条件分布——对许多复杂分布不可行;(2)坐标轴方向探索——对强相关变量效率低(见例3.7.4);(3)不适合离散+连续混合空间。\textbf{选择原则}:若条件分布简单(如共轭先验),优选Gibbs;否则用MH或两者混合。}

{\color{red}\textbf{[一次更新一个坐标的几何图像]:} Gibbs沿坐标轴逐一更新——想象在二维平面上,每步只能水平或竖直移动,禁止对角线!这种``轴对齐''限制有利有弊:优点是每个方向的移动由条件分布精确控制;缺点是对旋转的椭圆形分布(强相关变量),需要大量小步才能沿主轴移动,效率极低。对比MH可任意方向移动,更灵活但需接受/拒绝修正。}


设状态用 \(\mathbf{x} = \left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)\) 表示,设目标分布为 \(\pi \left( \mathbf{x}\right)\) ,用 \({\mathbf{x}}_{\left( -i\right) }\) 表示 \(\left( {{x}_{1},\ldots ,{x}_{i - 1}}\right.\) , \(\left. {{x}_{i + 1},\ldots ,{x}_{n}}\right)\) ,假设 \(\pi \left( \cdot \right)\) 的条件分布 \(p\left( {{x}_{i} \mid  {\mathbf{x}}_{\left( -i\right) }}\right)\) 都能够比较容易地抽样。

Gibbs 抽样每一步从条件分布中抽样, 可以轮流从每一分量抽样, 这样的算法称为系统扫描 Gibbs 抽样算法:

从 \(\pi \left( \mathbf{x}\right)\) 的取值区域中任意取一个初值 \({\mathbf{X}}^{\left( 0\right) }\)

for \(\left( {t\text{ in }0 : \left( {N - 1}\right) }\right) \{\)

for \(\left( {i\text{ in }1 : n}\right) \{\)

从条件分布 \(p\left( {{x}_{i} \mid  {X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)\) 中抽取 \({X}_{i}^{ * }\)

\}

令 \({\mathbf{X}}^{\left( t + 1\right) } \leftarrow  \left( {{X}_{1}^{ * },\ldots ,{X}_{n}^{ * }}\right)\)

\}

从条件分布抽样的次序也可以是随机选取各个分量,这样的算法称为随机扫描 Gibbs 抽样算法:

从 \(\pi \left( \mathbf{x}\right)\) 的取值区域中任意取一个初值 \({\mathbf{X}}^{\left( 0\right) }\)

for \(\left( {t\text{ in }0 : \left( {N - 1}\right) }\right) \{\)

按概率 \(\mathbf{\alpha } = \left( {{\alpha }_{1},\ldots ,{\alpha }_{n}}\right)\) 从 \(\left( {1,\ldots ,n}\right)\) 中随机抽取下标 \(i\)

从条件分布 \(p\left( {{x}_{i} \mid  {\mathbf{X}}_{\left( -i\right) }^{\left( t\right) }}\right)\) 中抽取 \({X}_{i}^{ * }\) ,令 \({\mathbf{X}}^{\left( t + 1\right) } = \left( {{X}_{1}^{\left( t\right) },\ldots ,{X}_{i - 1}^{\left( t\right) },{X}_{i}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)\)

\}

其中下标的抽样概率 \(\alpha\) 为事先给定。

容易看出,无论采用系统扫描还是随机扫描的 Gibbs 抽样,如果 \({\mathbf{X}}^{\left( t\right) }\) 服从目标分布,则 \({\mathbf{X}}^{\left( t + 1\right) }\) 也服从目标分布。以系统扫描方法为例,设在第 \(t + 1\) 步已经抽取了 \({X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * }\) ,令 \(\mathbf{Y} = \left( {{X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)\) ,设 \(\mathbf{Y} \sim  \pi \left( \cdot \right)\) 。下一步从 \(\pi \left( \cdot \right)\) 的边缘密度 \(p\left( {{x}_{i} \mid  {X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)\) 抽取 \({X}_{i}^{ * }\) ,则 \({\mathbf{Y}}^{ * }\overset{\bigtriangleup }{ = }\left( {{X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)\) 的分布密度在 \({\mathbf{Y}}^{ * }\) 处的值为

\[
p\left( {{X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)
\]

\[
= p\left( {{X}_{i}^{ * } \mid  {X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right) p\left( {{X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)
\]

\[
= \pi \left( {{X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)
\]

即 \(\mathbf{Y} \sim  \pi \left( \cdot \right)\) 则 \({\mathbf{Y}}^{ * } \sim  \pi \left( \cdot \right)\) 。

例 3.7.4. 设目标分布为二元正态分布,设 \(\mathbf{X} \sim  \pi \left( \mathbf{x}\right)\) 为

\[
\mathrm{N}\left\{  {\left( \begin{array}{l} 0 \\  0 \end{array}\right) ,\left( \begin{array}{ll} 1 & \rho \\  \rho & 1 \end{array}\right) }\right\}
\]

采用系统扫描 Gibbs 抽样方案, 每一步的迭代为,

\[
\text{ 抽取 }{X}_{1}^{\left( t + 1\right) } \mid  {X}_{2}^{\left( t\right) } \sim  \mathrm{N}\left( {\rho {X}_{2}^{\left( t\right) },1 - {\rho }^{2}}\right)
\]

\[
\text{ 抽取 }{X}_{2}^{\left( t + 1\right) } \mid  {X}_{1}^{\left( t + 1\right) } \sim  \mathrm{N}\left( {\rho {X}_{1}^{\left( t + 1\right) },1 - {\rho }^{2}}\right)
\]

递推可得

\[
\left( \begin{array}{l} {X}_{1}^{\left( t\right) } \\  {X}_{2}^{\left( t\right) } \end{array}\right)  \sim  \mathrm{N}\left\{  {\left( \begin{matrix} {\rho }^{{2t} - 1}{X}_{2}^{\left( 0\right) } \\  {\rho }^{2t}{X}_{2}^{\left( 0\right) } \end{matrix}\right) ,\left( \begin{matrix} 1 - {\rho }^{{4t} - 2} & \rho  - {\rho }^{{4t} - 1} \\  \rho  - {\rho }^{{4t} - 1} & 1 - {\rho }^{4t} \end{matrix}\right) }\right\}   \tag{3.96}
\]

当 \(t \rightarrow  \infty\) 时, \(\left( {{X}_{1}^{\left( t\right) },{X}_{2}^{\left( t\right) }}\right)\) 的期望与目标分布期望之差为 \(O\left( {\left| \rho \right| }^{2t}\right)\) ,方差与目标分布方差之差为 \(O\left( {\left| \rho \right| }^{4t}\right)\) 。

{\color{red}\textbf{[收敛速度的关键因素]:} 误差按$|\rho|^{2t}$指数衰减——相关系数$\rho$越接近1收敛越慢!当$\rho=0.99$时需要约460步才能使误差减小到1\%。这揭示了MCMC的``软肋'':强相关的变量导致缓慢混合。}

{\color{red}\textbf{[为什么相关性影响收敛]:} 物理解释:高相关($\rho \approx 1$)意味着目标分布是``细长''的椭圆形,主轴严重倾斜。Gibbs只能沿坐标轴移动(水平或竖直),要沿主轴方向前进,必须走``之字形''路径——水平一小步、竖直一小步,反复横跳!每步位移的主轴方向分量微小,故需要大量步数。对比:若$\rho=0$(圆形分布),坐标轴与主轴对齐,Gibbs高效。\textbf{解决方案}:参数重整化,将变量变换为近似独立的坐标系。}

{\color{red}\textbf{[二元正态的显式解]:} 例3.7.4的精确解(3.96)是Gibbs收敛性的罕见显式公式!从中可见:(1)均值以$O(\rho^{2t})$收敛到0;(2)方差以$O(\rho^{4t})$收敛到单位矩阵。初值$(0, X_2^{(0)})$的影响呈指数衰减。这个例子虽简单,但揭示了一般MCMC收敛的几何级数律——谱间隙(spectral gap,$1-|\lambda_2|$,其中$\lambda_2$是第二大特征值)决定收敛速度。$\rho$大→谱间隙小→收敛慢。}

例 3.7.5. 设目标分布为

\[
\pi \left( {x,y}\right)  \propto  \left( \begin{array}{l} n \\  x \end{array}\right) {y}^{x + \alpha  - 1}{\left( 1 - y\right) }^{n - x + \beta  - 1},x = 0,1,\ldots ,n,0 \leq  y \leq  1, \tag{3.97}
\]

则 \(X\left| {Y \sim  \mathrm{B}\left( {n,y}\right) ,Y}\right| X \sim  \operatorname{Beta}\left( {x + \alpha ,n - x + \beta }\right)\) 。易见 \(Y\) 的边缘分布为 \(\operatorname{Beta}\left( {\alpha ,\beta }\right)\) 。可以用 Gibbs 抽样方法模拟生成(X, Y)的样本链。

例 3.7.6. 在 Gibbs 抽样中, 每次变化的可以不是单个的分量, 而是两个或多个分量。例如, 设某个试验有 \(r\) 种不同结果,相应概率为 \(\mathbf{p} = \left( {{p}_{1},\ldots ,{p}_{r}}\right)\) (其中 \(\left. {\mathop{\sum }\limits_{{i = 1}}^{r}{p}_{i} = 1}\right)\) ,独立重复试验 \(n\) 次,各个结果出现的次数 \(\mathbf{X} = \left( {{X}_{1},\ldots ,{X}_{r}}\right)\) 服从多项分布。设 \(A = \left\{  {{X}_{1} \geq  1,\ldots ,{X}_{r} \geq  1}\right\}\) , 假设 \(P\left( A\right)\) 概率很小,要在条件 \(A\) 下对 \(\mathbf{X}\) 抽样,如果先生成 \(\mathbf{X}\) 的无条件样本再舍弃不符合条件 \(A\) 的部分则效率太低,可以采用如下的 Gibbs 抽样方法。

首先,任取初值 \({\mathbf{X}}^{\left( 0\right) }\) ,如 \({\mathbf{X}}^{\left( 0\right) } = \left( {1,\ldots ,1}\right)\) 。假设已生成了 \({\mathbf{X}}^{\left( t\right) }\) ,下一步首先从 \(\left( {1,\ldots ,r}\right)\) 中随机地抽取两个下标(i, j),令 \(s = {X}_{i}^{\left( t\right) } + {X}_{j}^{\left( t\right) }\) ,在给定 \({X}_{k},k \neq  i,j\) 为 \({\mathbf{X}}^{\left( t\right) }\) 对应元素的条件下, \(\left( {{X}_{i},{X}_{j}}\right)\) 的条件分布实际是 \(\left( {{X}_{i},{X}_{j}}\right)\) 在 \({X}_{i} + {X}_{j} = s\) 以及 \({X}_{i} \geq  1,{X}_{j} \geq  1\) 条件下的分布。于是,在以上条件下, \({X}_{i}\) 服从 \(\mathrm{B}\left( {s,{p}_{i}/\left( {{p}_{i} + {p}_{j}}\right) }\right)\) 分布限制在 \(1 \leq  {X}_{i} \leq  s - 1\) 条件下的分布,即

\[
{q}_{k}\overset{\bigtriangleup }{ = }P\left( {{X}_{i} = k \mid  {X}_{i} + {X}_{j} = s,{X}_{i} \geq  1,{X}_{j} \geq  1}\right)
\]

\[
= \frac{{C}_{s}^{k}{\left( \frac{{p}_{i}}{{p}_{i} + {p}_{j}}\right) }^{k}{\left( \frac{{p}_{j}}{{p}_{i} + {p}_{j}}\right) }^{s - k}}{1 - {\left( \frac{{p}_{j}}{{p}_{i} + {p}_{j}}\right) }^{s} - {\left( \frac{{p}_{i}}{{p}_{i} + {p}_{j}}\right) }^{s}},k = 1,2,\ldots ,s - 1.
\]

要生成这样的 \({X}_{i}\) 的抽样只要用生成离散型随机数的逆变换法。设抽取的 \({X}_{i}\) 值为 \({X}_{i}^{ * }\) ,取 \(\left( {{X}_{i}^{\left( t + 1\right) },{X}_{j}^{\left( t + 1\right) }}\right)  = \left( {{X}_{i}^{ * },s - {X}_{i}^{ * }}\right)\) ,取 \({\mathbf{X}}^{\left( t + 1\right) }\) 的其它元素为 \({\mathbf{X}}^{\left( t\right) }\) 的对应元素。如此重复就可以生成所需的 \(\mathbf{X}\) 在条件 \(A\) 下的抽样链。

\subsection{MCMC 计算软件 *}

{\color{red}\textbf{[MCMC的常见陷阱]:} (1)\textbf{虚假收敛}:轨迹图看似平稳,实则困在局部模式,遗漏主要质量。对策:多链从极端初值开始,观察是否收敛到同一区域;(2)\textbf{长程相关}:样本看似独立,实则存在长程相关(slow mixing),有效样本量远小于链长。对策:检查自相关图,确保在合理滞后(如lag 50)衰减到0;(3)\textbf{参数不可识别}:后验在某方向``平坦无界'',链无限漫游。对策:检查参数轨迹是否无明确中心,考虑重参数化或增强先验;(4)\textbf{过早停止}:误以为已收敛,实则仍在老化期。对策:丢弃至少前50\%样本,用多种诊断确认。}

{\color{red}\textbf{[何时信任MCMC结果]:} 满足以下所有条件时,可以合理信任:(1)多条链(≥3)从过度分散初值出发,都收敛到相同分布($\hat{R} < 1.1$);(2)轨迹图呈现稳定的``毛毛虫''状,无趋势或周期;(3)自相关在合理滞后内衰减,有效样本量足够大($n_{\text{eff}} > 1000$);(4)后验对先验不过度敏感(敏感性分析);(5)对真实数据,后验预测与观测数据一致(后验预测检验)。缺一不可!MCMC不是魔法,垃圾进垃圾出。}

{\color{red}\textbf{[MCMC计算成本]:} 典型应用的计算量:(1)\textbf{简单模型}(10维参数,简单似然):数千步,秒级;(2)\textbf{中等模型}(100维,需数值求解ODE):数万步,分钟到小时级;(3)\textbf{复杂模型}(高维层次模型,大数据):数十万步,小时到天级。瓶颈通常是似然函数计算,而非MCMC算法本身。优化策略:向量化计算、并行链、GPU加速、近似似然(如变分贝叶斯)。权衡:精确性vs计算成本。}

MCMC 是贝叶斯统计计算中最常用的计算工具。OpenBUGS 是一个成熟的 MCMC 计算开源软件 (见 Lunn et al(2009) \({}^{\left\lbrack  {29}\right\rbrack  }\) , Cowles(2013) \({}^{\left\lbrack  {16}\right\rbrack  }\) ,另一个类似的有关软件是 WinBUGS \({}^{\left\lbrack  {30}\right\rbrack  }\) ), 能够进行十分复杂的贝叶斯模型的计算,可以在 \(\mathrm{R}\) 中直接调用 OpenBUGS 进行计算。

{\color{red}\textbf{[实际应用]:} MCMC彻底改变了贝叶斯统计——1990年代前,复杂后验分布几乎无法计算;现在,成百上千维的层次模型都能轻松处理。从基因组学到气候模型,从经济预测到药物试验,MCMC无处不在。}

{\color{red}\textbf{[现代MCMC工具生态]:} OpenBUGS/WinBUGS虽经典但已过时。\textbf{现代工具}:(1)\textbf{Stan}(基于Hamiltonian MC)——自动求导,适合连续参数,诊断完善,推荐首选;(2)\textbf{PyMC}(Python)——灵活,与NumPy/JAX集成,社区活跃;(3)\textbf{JAGS}——类似BUGS但更现代,支持更多分布;(4)\textbf{TensorFlow Probability/Pyro}——与深度学习框架结合,适合大规模问题。选择原则:连续模型用Stan,离散/混合模型用PyMC/JAGS,大数据用TFP。}

OpenBUGS 采用 Gibbs 抽样方法从贝叶斯后验分布中抽样, 用户只需要指定先验分布和似然函数以及观测数据、已知参数, 以及并行地生成多少个马氏链、链的一些初值、运行步数。软件自动计算 Gibbs 抽样所需的条件分布, 产生马氏链, 并可以用图形和数值辅助判断收敛性, 给出后验推断的概括统计。

在 \(\mathrm{R}\) 中通过 BRugs 包调用 OpenBUGS 的功能。BRugs 用三类输入文件指定一个贝叶斯模型, 第一类文件指定似然函数和参数先验密度, 第二类文件指定已知参数、样本值, 第三类文件指定马氏链初值(并行产生多个链时需要指定多组初值)。 \(\mathrm{R}\) 的 coda 软件包可以帮助对 MCMC 抽样结果进行分析和诊断。下面用一个简单例子介绍在 \(\mathrm{R}\) 中用 BRugs 和 OpenBUGS 从贝叶斯后验中抽样的基本步骤。

\iffalse
例 3.7.7. 对例3.7.2,用 \(\mathrm{R}\) 的 BRugs 包调用 OpenBUGS 来计算。设 \(\left( {{X}_{1},\ldots ,{X}_{5}}\right)\) 的观测值为(74,85,69,17,5)。

OpenBUGS 用模型文件描述随机变量分布和对参数的依赖关系, 以及参数之间的关系。 首先建立如下模型文件,保存在文件 pfl-model.txt 中:

\HRule

model

\{

\hspace*{1em} \(p\left\lbrack  1\right\rbrack   <  - 1/3\)

\hspace*{1em} \(p\left\lbrack  2\right\rbrack   <  - \left( {1 - b}\right) /3\)

\hspace*{1em} \(p\left\lbrack  3\right\rbrack   <  - \left( {1 - 2 * b}\right) /3\)

\hspace*{1em} \(p\left\lbrack  4\right\rbrack   <  - 2 * b/3\)

\hspace*{1em} \(p\left\lbrack  5\right\rbrack   <  - b/3\)

\hspace*{1em} b <- b2 / 2

\hspace*{1em} b2 \~ dbeta(1,1)

\hspace*{1em} x[1:5] \~ dmulti(p[1:5], N)

\}

\HRule

文件中用向左的箭头表示确定性的关系,用方括号加序号表示下标,用 \(\sim\) 表示左边的变量服从右边的分布。这里, \(b\) 为参数 \(\beta ,\mathrm{b}2 = {2b}\) 服从 \(\operatorname{Beta}\left( {1,1}\right)\) 分布即 \(\mathrm{U}\left( {0,1}\right)\) 分布,于是 \(\beta\) 有先验分布 \(\mathrm{U}\left( {0,{0.5}}\right)\) 。 \(\mathrm{x}\left\lbrack  {1 : 5}\right\rbrack\) 表示向量 \(\left( {{x}_{1},\ldots ,{x}_{5}}\right)\) ,服从多项分布,参数为 \(\left( {{p}_{1},\ldots ,{p}_{5}}\right)\) ,试验次数为 \(N\) 。 \(N\) 在模型文件中没有指定,将在数据文件中给出。

建立如下的数据文件, 保存在文件 pfl-data.txt 中:

\HRule

list(x=c(74, 85, 69, 17, 5),

\hspace*{1em} N=250)

\HRule

这样的数据文件是 \(\mathrm{R}\) 软件的列表格式,列表中的标量和向量为通常的 \(\mathrm{R}\) 程序写法,矩阵用如

\HRule

list(M=structure(.Data=c(1,2,3,4,5,6),.Dim=c(3,2)))

\HRule

表示,其中. Dim 给出矩阵的行、列数 \(\left( {3 \times  2}\right)\) ,. Data 给出按行排列的所有元素(第一行为 1,2,第二行为3,4,第三行为5,6)。

OpenBUGS 需要用户指定各个链的要抽样的参数的初值, 这里我们要抽样的参数是 b, 但 \(\mathrm{b}\) 是由 \(\mathrm{b}2\) 计算得到的,所以对 \(\mathrm{b}2\) 设置初值。设有如下两个初值文件,不同链的初值应

尽可能不同, 文件 pfl-inits1.txt 内容为:

\HRule

list(b2=0.2)

\HRule

文件 pfl-inits2.txt 内容为:

\HRule

list(b2=0.8)

\HRule

在 \(\mathrm{R}\) 中,首先调用 BRugs 软件包:

\HRule

require (BRugs)

\HRule

然后, 读入并检查模型文件:

\HRule

modelCheck('pfl-model.txt')

\HRule

读入数据文件:

\HRule

modelData ( ' pfl-data.txt ' )

\HRule

下面, 对模型和数据进行编译, 得到抽样方案, 下面的语句要求并行运行两个链:

\HRule

modelCompile(numChains=2)

\HRule

准备迭代地生成 MCMC 抽样了, 首先指定初值:

\HRule

modelInits(c('pfl-inits1.txt', 'pfl-inits2.txt'))

\HRule

下面, 先试验性地运行 1000 次, 作为老化期:

\HRule

modelUpdate (1000)

\HRule

现在才指定抽样要输出那些随机变量的随机数:

\HRule

samplesSet \(\left( {c\left( {{}^{\prime }{b}^{\prime },{}^{\prime }p\left\lbrack  {1 : 5}\right\rbrack  {}^{\prime }}\right) }\right)\)

\HRule

现在可以抽样了, 指定运行 10000 次, 两个链并行运行:

\HRule

modelUpdate (10000)

\HRule

得到抽样后,可以把抽样的结果保存在 \(\mathrm{R}\) 的变量中,比如

\HRule

b2chains <- samplesHistory ( 'b', plot=FALSE )

\HRule

得到一个列表,有唯一的元素 \(\mathrm{b}\) ,为 \(2 \times  {10000}\) 的矩阵,每行是一个链的记录。sampleHistory 可以抽样链的曲线图, 如

samplesHistory('b')

OpenBUGS 提供了一系列的简单统计和收敛诊断功能。如下程序列出各抽样变量的简单统计:

\HRule

print ( samplesStats ( ``*'' ) )

\HRule

结果为:

\HRule

\hspace*{1em} <table><tr><td/><td>mean</td><td>sd</td><td>MC\_error</td><td>val2.5pc</td><td>median</td><td>val97.5pc</td><td>start</td><td>sample</td></tr><tr><td>b</td><td>0.08757</td><td>0.016900</td><td>1.168e-04</td><td>0.05737</td><td>0.08668</td><td>0.12330</td><td>1001</td><td>20000</td></tr><tr><td> \(\mathrm{p}\left\lbrack  2\right\rbrack\) </td><td>0.30410</td><td>0.005632</td><td>3.892e-05</td><td>0.29230</td><td>0.30440</td><td>0.31420</td><td>1001</td><td>20000</td></tr><tr><td> \(\mathrm{p}\left\lbrack  3\right\rbrack\) </td><td>0.27500</td><td>0.011260</td><td>7.784e-05</td><td>0.25120</td><td>0.27550</td><td>0.29510</td><td>1001</td><td>20000</td></tr><tr><td>p [4]</td><td>0.05838</td><td>0.011260</td><td>7.784e-05</td><td>0.03824</td><td>0.05778</td><td>0.08217</td><td>1001</td><td>20000</td></tr><tr><td> \(\mathrm{p}\left\lbrack  5\right\rbrack\) </td><td>0.02919</td><td>0.005632</td><td>3.892e-05</td><td>0.01912</td><td>0.02889</td><td>0.04109</td><td>1001</td><td>20000</td></tr></table>

\HRule

统计使用所有链的数据。可以看出, \(\beta\) 的后验均值为 0.08757。表中的 \({\mathrm{{MC}}}_{ - }\) error 表示估计后验均值时由于随机模拟导致的误差的标准差的估计, 这个标准差估计针对抽样自相关性进行了校正。在老化期之后的运行次数越多 MC\_error 越小, 一个常用的经验规则是保证 MC\_error 小于后验标准差的 \(5\%\) ,这里的结果提示还需要更多的运行次数。

{\color{red}\textbf{[蒙特卡洛误差]:} MC\_error与普通标准误的区别在于``相关性校正''——独立样本的标准误是$\sigma/\sqrt{n}$,但MCMC样本相关,有效样本量远小于$n$。如果自相关强,可能10000个样本的信息量只等价于100个独立样本!}

val \({2.5}\mathrm{{pc}}\) 和 val \({97.5}\mathrm{{pc}}\) 是抽样的后验分布的 2.5\% 和 97.5\% 分位数的估计值,由此得到 \(\beta\) 的水平 95\% 的可信区间 (credible interval) 为(0.05734,0.1233)。

如下程序画出 \(\mathrm{b}\) 的后验密度估计:

\HRule

samplesDensity ( ' b ', mfrow=c ( 1 , 1 ) )

\HRule

如下程序画抽样的 \(\mathrm{b}\) 的自相关函数估计:

\HRule

samplesAutoC('b', 1, mfrow=c(1,1))

\HRule

当自相关函数很大而且衰减缓慢时生成的抽样链的效率较低。本例中的自相关函数基本表现为不相关列。

为了检查链是否收敛, OpenBUGS 提供了 BGR 统计量图:

\HRule

samplesBgr ( ' b ', mfrow=c (1,1 ) )

\HRule
\fi
BGR 统计量的原理是考虑并行运行的每个链内部的变化情况, 以及把所有的链合并在一起的变化情况, 对这两种变化情况进行比较, 当链收敛时, 每个链内部的变化情况应该与合并在一起的变化情况很类似。类似于单因素方差分析中组间平方和与组内平方和的比较。当 BGR 图中的红色线接近于 1 并且三条线都保持稳定时就提示链收敛了。

{\color{red}\textbf{[收敛诊断的哲学]:} BGR(Gelman-Rubin)统计量体现``多链共识''——如果不同起点的链都到达同一分布,说明确实收敛了;如果链间差异大于链内差异,说明还在``各走各路''。这类似政治民主:多方达成共识才可靠!}

{\color{red}\textbf{[BGR统计量的数学定义]:} 设有$m$条链,每条长度$n$(丢弃老化期后)。定义:(1)\textbf{链内方差}$W = \frac{1}{m}\sum_{j=1}^m s_j^2$,其中$s_j^2$是第$j$条链的样本方差——衡量单链内涨落;(2)\textbf{链间方差}$B = \frac{n}{m-1}\sum_{j=1}^m (\bar{\theta}_j - \bar{\theta})^2$,其中$\bar{\theta}_j$是第$j$链均值,$\bar{\theta}$是总均值——衡量链间差异;(3)\textbf{$\hat{R}$统计量}$= \sqrt{\frac{\hat{V}}{W}}$,其中$\hat{V} = \frac{n-1}{n}W + \frac{1}{n}B$是总方差估计。\textbf{判据}:$\hat{R} \lesssim 1.1$认为收敛(越接近1越好)。}

{\color{red}\textbf{[BGR为何有效]:} 直觉:(1)收敛前:链还在``爬山'',从不同初值出发的链位于不同位置→链间方差$B$大,远超链内方差$W$→$\hat{R} \gg 1$;(2)收敛后:所有链都在平稳分布$\pi$附近涨落,位置相似→$B \approx W$→$\hat{R} \approx 1$。类比ANOVA:若组间差异显著大于组内差异,说明``组''这个因素有效(这里意味着初值仍有影响,未收敛);若组间组内差不多,说明``组''无关紧要(初值被遗忘,已收敛)。}

{\color{red}\textbf{[收敛诊断的实践建议]:} (1)\textbf{多链策略}:至少跑3-5条链,从\textbf{过度分散}的初值开始(覆盖参数空间的不同区域);(2)\textbf{多种诊断}:不要只看BGR!结合轨迹图(trace plot)、自相关图、有效样本量等;(3)\textbf{保守原则}:收敛诊断只能识别``明显未收敛'',无法证明``确实收敛''——可能存在链未发现的远距离多峰!谨慎对待临界$\hat{R} \approx 1.1$的情况;(4)\textbf{参数异质性}:不同参数收敛速度不同,需对\textbf{所有}关心的参数检查,最慢的决定总老化期。}
